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Abstract — A variational approach for filling-in regions of 
missing data in digital images is introduced in this paper. The ap- 
proach is based on joint interpolation of the image gray-levels and 
gradient/isophotes directions, smoothly extending in an automatic 
fashion the isophote lines into the holes of missing data. This in- 
terpolation is computed by solving the variational problem via its 
gradient descent flow, which leads to a set of coupled second order 
partial differential equations, one for the gray-levels and one for 
the gradient orientations. The process underlying this approach 
can be considered as an interpretation of the Gestaltist* s principle 
of good continuation. No limitations are imposed on the topology 
of the holes, and all regions of missing data can be simultaneously 
processed, even if they are surrounded by completely different 
structures. Applications of this technique include the restoration 
of old photographs and removal of superimposed text like dates, 
subtitles, or publicity. Examples of these applications are given. 
We conclude the paper with a number of theoretical results on 
the proposed variational approach and its corresponding gradient 
descent flow. 

Index Terms— Filling-in, Gestalt principles, image gradients, 
image gray-levels, interpolation, partial differential equations, 
variational approach. 

I. Introduction 

FILLING-IN missing data in digital images has a number of 
fundamental applications. They range from removing ob- 
jects from a scene all the way to retouching damaged paintings 
and photographs. The basic idea is to fill-in the gap of missing 
data in a form that it is nondetectable by an ordinary observer 
In art, this process is called inpainting [6], [7], [16], [25], [39]. 

Since the early days of art and photography, filling-in and in- 
painting has been done by professional artists. Imitating their 
performance with semi-automatic digital techniques is currently 
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an active area of research; see, for example, [8], [23], [27], 
and [28], and the papers discussed in this paper. The goal of 
this work is to introduce a novel algorithm for automatically 
filling-in gaps in the image. In this paper, we follow the sugges- 
tions in the conclusions section in [6] and introduce an energy 
functional based on an interpretation of Gestaltist's good contin- 
uation principle. Suppose that we are given an image u 0 : R — > 
[a.b], where R is a square of JR 2 and a < b, and 12 is an open 
bounded subset of R with Lipschitz continuous boundary. We 
shall call ft the hole or gap. We want to fill-in the hole O based 
on the geometric and photometric information outside the hole. 
For that we use what we call a band around Q, i.e., we consider 
an open region ft of R such that Q C Q (ft is the closure of the 
set). The band we refer to will be the set B — ft\ft. To fill-in 
the hole ft we use the information of contained in B, mainly 
gray level and the vector field of isophotes (level sets) directions 
of uo in B. 1 We attempt to continue the level sets of u in B in- 
side ft taking into account the principle of good continuation. 
We propose an energy functional which takes into account these 
principles interpreted in a suitable way. The energy functional 
we propose has to be minimized with respect to two variables: 
a vector field 0 which represents the directions of the level lines 
of u, and the gray level u. 0 and u are constrained in the band 
B by their known values there. The use of the vector field of 
directions 0 is one of the main points of the algorithm presented 
in this paper, which permits the level sets to smoothly continue 
inside the hole. We are then continuing both the geometric and 
photometric properties of the image inside the hole. 

Let us finally say that the only user interaction required by the 
algorithm here introduced is to mark the regions to be filled-in. 
Although a number of techniques exist for the semi-automatic 
detection of image defects (mainly in films), addressing this is 
not part of the scope of this paper. Since the algorithm here 
presented can be used not just to restore damaged photographs 
but also to remove undesired objects and writings on the image, 
the regions must be marked by the user, since they depend on 
his/her subjective selection. 

A. Closely Related Approaches 

Before proceeding with the detailed description of our algo- 
rithm, let us comment on related work. Note that image de- 
noising is different to filling-in, since the regions of missing 
data are usually large. That is, regions occupied by top to bottom 
scratches along several film frames, long cracks in photographs, 

'The width of the band is such that it conveys the boundary information, 
mainly the gray level and isophotes direction, and numerically it depends on 
the minimal stencil needed by the implementation. 
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superimposed large fonts, and so on, are of significant larger size 
than the type of noise assumed in common image enhancement 
algorithms. In addition, in common image enhancement appli- 
cations, the pixels contain both information about the real data 
and the noise (e.g., image plus noise for additive noise), while in 
our application there is no significant information in the region 
to be inpainted. 

A very active area related to the work here presented is the 
restoration of damaged films. The basic idea here is to use in- 
formation from past and future frames to restore the current one, 
e.g., [23], [28]. Of course, this general approach cannot be used 
when dealing with still images. 

Another area related to the work here described is texture 
synthesis. The basic idea here is to select a texture and synthe- 
size it inside the region to be filled-in (the hole). Although out- 
standing texture synthesis results have been reported in the lit- 
erature, e.g., [14], [19], [22], and [36], these algorithms require 
the user to select the texture to be copied into the hole. For im- 
ages in which the region to be replaced covers several different 
structures, the user would need to go through the tremendous 
work of segmenting them and searching corresponding replace- 
ments throughout the picture. Although part of this search can be 
performed automatically, the overall process is extremely time 
consuming and requires the nontrivial selection of many critical 
parameters, e.g., [14]. 

Last, a number of fundamental works on the topic of disocclu- 
sion and line continuation have been reported in the literature, 
and these are the closest to our approach. A pioneering contri- 
bution in this area is described in [33]. The authors presented a 
technique for removing occlusions with the goal of image seg- 
mentation. Since the region to be filled-in can be considered as 
occluding objects, removing occlusions is analogous to image 
inpainting. The basic idea suggested by the authors is to con- 
nect T-junctions at the occluding boundaries of objects with 
elastica minimizing curves (see later in this paper for the exact 
definition of elastica curves). The technique was primarily de- 
veloped for simple images obtained from a segmentation, with 
only a few objects with constant gray-levels. Thus, they ended 
up by connecting T-junctions at the same gray level. (Other re- 
searchers, e.g., Jacobs et al. have followed this interesting re- 
search area, mainly developing techniques for smooth curve 
continuation.) Masnou and Morel [31], [32] recently extended 
these ideas, presenting a very elegant and inspiring formal vari- 
ational formulation for disocclusion and a particular practical 
algorithm implementing some of the ideas in this formulation. 
The algorithm fills- in by joining with geodesic curves the points 
of the isophotes arriving at the boundary of the region to be 
inpainted. The holes in their algorithm are limited to having 
simple topology. In addition, the angle with which the level lines 
arrive at the boundary of the holes are not (well) preserved, and 
the algorithm uses straight lines to join equal gray value pixels. 
The present work was motivated, in part, by their work. In some 
sense, as shown by the work of Masnou and Morel [3 1 ], [32] and 
this paper, operators as the one proposed in [33] can be used to 
first interpolate level-sets, and then reconstruct from them the 
interpolated gray- level image. 

Recently, we have addressed the concept of smooth contin- 
uation of information in the level-lines direction in [6], [7]. 



We proposed an algorithm, inspired by partial differential equa- 
tions, that propagates the image Laplacian in this direction. The 
algorithm attempts to imitate basic approaches used by pro- 
fessional restorators. The algorithm also introduces the impor- 
tance of propagating both the gradient direction (geometry) and 
gray-values (photometry) of the image in a band surrounding 
the hole to be filled-in. It is part of the goal of the current paper 
to adopt some of the ideas of [6] and [7], while deviating from 
the particular model in order to be able to define a formal vari- 
ational approach to the filling-in problem. 

The work in [6] and [7] inspired a very elegant approach to the 
filling-in problem recently reported in [9] (this work was per- 
formed independently to the one reported in this paper). 2 The 
authors present a clear and intuitive axiomatic approach to the 
problem. The main algorithm they propose after an interesting 
discussion of the inpainting problem is to minimize the total 
variation (TV), [35], of the image inside the hole (they also use, 
as proposed in [6], [7] and further studied in this paper, a band 
surrounding the region). They address in addition the interpola- 
tion and filling-in in the presence of noise, a very important ad- 
ditional contribution. As in the work of Masnou and Morel, their 
interpolation is limited to creating straight isophotes, not neces- 
sarily smoothly continued from the hole boundary, and mainly 
is developed (as the authors clearly state) for small holes. Al- 
though straight connections give visually pleasant results for 
small holes, it is important to develop a theory that permits inter- 
polation of level lines across large gaps, where connecting with 
straight lines will be unpleasant even for simple images. As we 
will argue later in this paper, in order to obtain such a smooth 
interpolation and continuation of isophotes, it is necessary to 
go into high-order partial differential equations (PDEs) or sys- 
tems of PDEs, as done in [6], [7] and here. Note that considering 
the angle of arrival of the level lines at the gap, and pursuing a 
smooth interpolation of it as done here, it is also supported by 
research in perception, from the Gestalt to more recent work, 
e.g., [34]. 

To conclude this section, the interested reader is referred to 
the works of Nitzberg-Mumford-Shiota, Masnou-Morel, and 
Chan-Shen (as well as our previous work) to study other inter- 
esting and very related techniques for filling- in. 

Let us conclude with the plan of the paper In Section II, 
we introduce the problem, the functional spaces and the energy 
functional for image inpainting. To clarify the meaning of the 
functional, we also discuss the particular case where we inter- 
polate the gray level, knowing the vector field of directions. Sec- 
tion III is devoted to numerical experiments. Section IV contains 
some conclusions. Finally, Appendices I and II are devoted, re- 
spectively, to the proof of existence of minimizers for the energy 
functional introduced in this paper and to details on the numer- 
ical implementation of the proposed equations. 

II. Joint Interpolation of Vector Fields 
and Gray Values 

Let uo R — > M be an image defined on a domain R of 
JR 2 , which we may suppose to be a square. Let O, f2 be two 

2 We thank the authors for providing us with a preliminary report of their work. 
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open bounded domains in JR 2 with Lipschitz boundary and sup- 
pose that ft C ft cc R. To simplify our presentation we shall 
assume that ft does not touch the boundary of the image do- 
main R. Let B := ft\ft. B will be called the band around ft. 
Suppose that a function uo is given in B which, for the mo- 
ment being, we shall assume to be smooth in the closure, B, of 
B (later we shall assume that no is of bounded variation, i.e., 
uq e BV(B)). Let 0 O be the vector field of directions of the 
gradient of w 0 on £?, i.e., Oq is a vector field with values in JR 2 
satisfying 6 0 {x) • Vua(x) = | Vw 0 (^)| and \0a(x)\ < 1 (ideally 
1 or 0, see below). This is the information we shall use on B. 

We pose the image inpainting problem in the following form: 
Can we extend (in a reasonable way) the pair of functions 
(uo. 9o) from the band ft\ft to a pair of functions (u, 0) defined 
inside ft? Of course, we will have to precise what we mean by 
a reasonable way. We shall discuss and analyze a variational 
formulation of this filling-in problem and discuss possible 
energy functionals, and their corresponding gradient descent 
flows, which give a solution to it. The data are given on the band 
B and we should constraint the solution (u, 9) to be near the 
data on B. The vector field 0 should satisfy |0| < 1 on ft and 
should be related to u by trying to impose that 9 • Vu — |Vu|, 
i.e., we should impose that 9 is related to the vector field of 
directions of the gradient of u. The condition \9(x)\ < 1 should 
be interpreted as a relaxation of this. Indeed, it may happen 
that 0(x) = 0 (flat regions) and then we cannot normalize 
the vector field to a unit vector. We should have in mind that 
the ideal case would be that 9 = Vu/|Vu|, u being a smooth 
function with Vu(x) / 0 for all x € ft. Finally, we should 
impose that the vector field 9o in the band tries to continue 
smoothly to 9 inside ft. We shall impose this by observing that, 
in case 9 represents the directions of the normals to the level 
lines of u, i.e., of the curves u(xi,x 2 ) — A, A 6 JR, then a term 
like div(0) represents its curvature. Motivated by the principle 
of smooth continuation, our energy functional should contain 
terms integrating div(0). Indeed, collecting all the observations 
above, we propose to minimize a functional of the form 

Minimize f\div(9)\ p {a + 6| Vk * u\) dx 
Jq 

+ a f(\Vu\-9-Vu)dx (1) 
Jq 

where 

a, 6, and a positive constants; 

A: smoothing kernel; 

0, u submitted to constraints. 

We need to give a sense to the integrals appearing in the above 
expression and to make precise the admissible class of functions 
where the functional has to be minimized. For that, we need to 
introduce some function spaces. This is done in the following 
section. Once this has been formally addressed, we will further 
discuss the underlying concepts of the above functional. 

A. Function Spaces 

Let us first recall the definition of B V functions and total vari- 
ation. Let Q be an open set. A function u € L l (Q) whose partial 



derivatives in the sense of distributions are measures with finite 
total variation in Q is called a function of bounded variation. 
The class of such functions will be denoted by BV(Q). Thus 
u e BV(Q) if and only if there are Radon measures pi , . . . , p N 
denned in Q with finite total mass in Q and 

/ uDi<pdx = - / ipdfa (2) 
Jq Jq 

for all <p e C^iQ), i = 1, . . . , N. Thus, the gradient of u is a 
vector valued measure with finite total variation 

|| V«|| =sup jj^ udivc^dx : <p e CS>{Q,lR n \ 

\y(x)\ < 1 forx e q| . (3) 

The space BV(Q) is endowed with the norm 

IMIbv = |HIlmq) + II v ^II- (4) 

We say that a measurable set E C Q has finite perimeter in 
Q if its indicator function \E € BV(Q). If u € BV(Q) al- 
most all its level sets [u > A] = {x e Q : u{x) > A} are 
sets of finite perimeter. For sets of finite perimeter E one can 
define the essential boundary d* E, which is rectifable with fi- 
nite H A -1 measure, and compute the normal to the level set at 
H N ~ l almost all points of d*E. Thus at almost all points of 
almost all level sets of u e BV(Q) we may define a normal 
vector 6{x). This vector field of normals 9 can be also denned 
(hence extended to all Q) as the Radon-Nikodym derivative of 
the measure Vu with respect to | Vu\, i.e., it formally satisfies 
0 - Via = |Vu| and, also, |0| < 1 a.e. For further information 
concerning functions of bounded variation we refer to [ 1 ], [ 1 7], 
and [40]. 

Let us now introduce the function spaces for 0. Let Q be an 
open bounded subset of JR 2 with a Lipschitz boundary. We de- 
fine 

W^idiv, Q) = {9e L p (Q) 2 : div(0) € L p (Q)}, 

1 < p < cc, 

and 

M(div, Q) = {6 e L X (Q) 2 : div(0) 

is a Radon measure in Q}. 

The Trace Theorem [2], [10] guarantees that the normal compo- 
nent 9'u\oq, is well defined for vector fields 9 in W^'f^div, Q), 
or in M (div, Q). To simplify our notation we shall assume that 
M /1 ' 1 (div, Q) represents the space M(div, Q). 

Next, we shall give a sense to the integrals of bounded vector 
fields with divergence in W integrated with respect to the gra- 
dient of a B V function. For that, we shall need some results from 
[2] (see also [11] and [26]). Let Q be an open bounded subset of 
M n with Lipschitz continuous boundary. Let p > 1 and q > 1 
be such that (1/p) + (1/q) = 1. Following [2], let 

X(Q) P = {ze L°°{Q, JR n ) : div(*) e If{Q)}. (5) 
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If z 6 X(Q) p and w e BV(Q) O L<*{Q) we define the func- 
tional 0, Vw) : qj°(9) — JR by the formula 

Vw), <p) = - / w(pdi\(z)dx- / wz-Vipdx. (6) 

Then, Vuy) is a Radon measure in Q 

/ (z,Vu;) = / z ■ Vwc& (7) 

for all iu€ W^QJni/^Q) and 

/ (z 5 vJ < / |(*,Vw)| < ||z||oo / ||V^|| (8) 
JB I Jb Jb 

for any Borel set B C Q. 

In [2], a weak trace on <9Q of the normal component of z € 
X(Q) P is defined. Concretely, it is proved that there exists a 
linear operator 7 : X(Q) P -> L°°(DQ) such that 

II7WIU < 
-y(z)(x) = z(z) • n(x) 

for all a? 6 if z € C 1 (Q, M N ) 

where n(x) denotes the outer unit normal at x € dQ. We shall 
denote f y(z)(x) by [z ) n](x). Moreover, the following Green's 
formula, relating the function [z,n] and the measure (z, Viu), 
for 2 e X{Q) P and w € BV(Q) n is established 

/ wdiv(*)dr + / (^Vw)= / [a? f n]iu<UT N " 1 . (9) 

If no confusion arises, we shall denote z • V w instead of (z y Vu) 
for z e -Y(Q) p ,u; e BV(Q) n These results will be 

used in the Proof of Theorem 2 given in Appendix I. 

B. Energy Derivation and Interpretation 

One of the key concepts above was the band around the hole. 
The band is of local character but in principle it could be ex- 
tended to all the known part of the image. Obviously, what hap- 
pens at distant parts can be independent or not from what hap- 
pens at the hole, but, in our construction below, we suppose that 
only a narrow band around the hole influences what happens in- 
side the hole. Could we fill-in without the band? To discuss this 
suppose that we are given the image of Fig. 1(a) which is a gray 
band on a black background partially occluded by a square i I. 
We suppose that the sides of the square hole ft are orthogonal 
to the level lines of the original image. In these conditions, the 
normal component of the vector field Oq outside ft is null at 5ft. 
Thus if the boundary data is just Oq * n \on, we would have that 
#o ■ n\on — 0- ln particular, the vector field 0 = 0 satisfies this 
condition. If we are not able to propagate 0 inside ft this may 
become an unpleasant situation, since this would mean that we 
do no propagate the values of u at the boundary. If we write the 
functional (1) with 0 — 0, it turns out to be the Total Variation 
[35]. The decision of extending the gray band or filling the hole 
with the black level would be taken as a function of the perimeter 
of the discontinuities of the function in the hole. Then the result 
of interpolating Fig. 1(a) using total variation would be that of 




(a) (b) (c) 



Fig. 1. (a) Original image, (b) filling-in via total variation-based techniques, 
and (c) filling-in via our proposed algorithm. 




Fig. 2. Region to be filled-in and corresponding surrounding band from where 
the information is considered. 

Fig. 1(b) and not the one in Fig. 1(c), because the interpolating 
lines in Fig. 1(b) are shorter than the ones in Fig. 1(c). To over- 
come this situation we introduce the band around the hole. The 
introduction of the band permits us to effectively incorporate in 
the functional the information given by the vector field 6 out- 
side ft. In Fig. 1 (b), we display the result of the interpolation 
with 9 = 0 on ft. In Fig. 1 (c), we show the result of the interpo- 
lation using the functional we shall completely describe below, 
which takes into account the band B and computes the vector 
field 0 in ft = ft U B. 

Thus, let B be a band around ft with a Lipschitz boundary 
containing the boundary of ft (see Fig. 2). As we made explicit 
above, B = ft\ft. Given the band B and the function u 0 of 
bounded variation in B, we define the space 

BV(ft,£,<u 0 ) = {u e BV(ft) : a = u 0 in B}. 

Let 0 O : B — > Hi 2 be a vector field of directions of the gradient 
of uo, i.e., 1 0o I < 1 and #0 * Vwo = V* u o\ as measures in B 
(therefore, a.e.). In practice we shall constraint the vector field 0 
to be the vector field of directions of u 0 only indirectly, through 
the functional. We could also introduce this as a constraint or 
with a penalty term J B (0 - 6q) 2 (see also [9], [35] for penalties 
of this form in the gray values). 

Combining the previous elements, the band, the relations be- 
tween 0 and u, and the smoothness term on 0, we propose to 
interpolate the pair (0, u) in ft by minimizing the functional 

Minimize / |div(0)| p (a + b\Vk * u\) dx 

+ a I (\Vu\-9-Vu)dx 
Jq 

6 e W l *{4iv,ti) 
u e BV(ft,£, u 0 ) 
0 . n \ d n = 9 0 • n|an 
1*1 <1 

M < ll«olL~(B) 0°) 
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where a, a > 0, b > 0, k denotes a regularizing kernel of class 
C 1 such that k(x) > 0 a.e. and n{x) denotes the outer unit 
normal at x G Oil. The previous functional is coercive and ad- 
mits a minimum in the class of functions described above if 
p > 1. The case p — 1 is under study. The functional can be 
interpreted as a formulation of the principle of good continua- 
tion and amodal completion as formulated in the Gestalt theory 
of vision. The following remarks contain heuristic arguments 
which may help to understand our choice. In next subsection 
we shall explain in more detail the role of the term coupling 6 
and u. 
Remarks: 

1) The constant b is >0. If u is the characteristic function of 
the region enclosed by a curve C then a term like 

[\div(9W\Vu\ (11) 

is related to J c \k\ p ds, where k is the Euclidean curvature 
(of the level-sets). If p = 2, this term appears in Eider's 
elastica 

J (a + (JK 2 )ds, a,/3>0. (12) 

Eider's elastica (12) was proposed in [33] as a technique 
for removing occlusions with the goal of image segmenta- 
tion, since this criterion yields smooth, short, and not too 
curvy curves. In terms of characteristic functions, Euler's 
elastica can be written as 

/|V«|(a + /i|div(^)|) 2 . (13) 

In [5], it was shown that this functional is not lower 
semicontinuous. The functional proposed by Masnou 
and Morel [3 1 ], [32] can be interpreted as a relaxation of 
it, since it integrates functional like the elastica (plus the 
angle that the curve makes with the corresponding level 
line arriving at the boundary) along the level lines of the 
function u. Our functional can be also considered as a 
relaxed formulation of the energy of the elastica. For that, 
we introduced 0 as a independent variable, and we tried 
to couple it to u by imposing that 0 * Vu = |Vw|. This 
restriction could be directly incorporated to the model as 
a constraint. We choose to incorporate this constraint as 
a penalization term 

f \Vu\ - 0 • V« (14) 
J o 

(see the next subsection for a detailed discussion of this 
term). Finally, let us say that for mathematical reasons we 
have convolved the Vu term of (1 1) to be able to prove 
the existence of a minimum for (10). From a theoretical 
point of view, this may invalidate our previous comments. 
But, from a practical point of view, it gives a weight to the 
curve of discontinuities of the image. 

2) The constant a has to be >0. Otherwise we do not get 
an LP bound on div(0). Now, let us comment on the two 



terms containing div(0). Heuristically, if we do not com- 
pute 0 in a proper way, in a continuous image like in 
Fig. 1, 0 could be zero except on a set of curves. Then 
0 = 0 a.e. on B (or on U) and a term like 

[\div{9)\ p dx (15) 

would produce a null value since div(0) — 0. On the other 
hand, a term like (11) would integrate a power of the cur- 
vature on the level line corresponding to the discontinuity 
of the image and it would guarantee that the functional 
is not null. This argument is only heuristic and not com- 
pletely justified. Indeed, we believe that in such example 
as in Fig. 1 , a term like (15) would induce a regularizing 
effect on 6 and the support of 0 would not be a curve any 
more. In that case, the integral (15) would not be null. 

3) Related to the question discussed in the last comment is 
the possibility to compute a regularized vector field of di- 
rections for images which are constant except at jump dis- 
continuities. A direct computation of the vector field out- 
side the hole in an image like Fig. 1(a) gives a null vector 
field at all points except the points on the level line sepa- 
rating the black from the white region. This may not be a 
good starting point to extend reasonably the vector field 
0 inside il To initialize the algorithm of steepest descent, 
a regularization of 9 outside Q, may be constructed as the 
vector field of directions of the image U 0 (t,x) obtained 
by regularizing ^0(^)7^ £ B> with 

du 

— = 0 mS = (0,oo) x 8B 
or) 

u(0,x) — uq(x) for x € B, (16) 

As it is shown in [3], this equation permits a regular- 
ization of the vector field of directions of the gradient 
of u, i.e., there is a vector field z.\z\ < 1, such that 
u t = div(z) and j Q z ■ Vu — J B |Vu|. Moreover, for 
each t > 0,div(*(t)) e D'(B) if u 0 6 W(B) for all 
p > 1. In this way, we initialize the steepest descent al- 
gorithm described in Section III with a regularized vector 
field 0. This again raises a question, namely, if this ad-hoc 
regularization is really needed or a regularization takes 
place with the algorithm itself, if we use an implicit nu- 
merical scheme to solve (10). 

4) The bound |it| < ||«o||l«(s) can be replaced by a con- 
stant depending on H'uoIIz^b) The constraint that u — 
uq in B could be relaxed by adding a penalty term like 
f B (u — tio) 2 Similarly, we could add a penalty term to 
constraint 0 to be near 0 O inside B. In this case, we should 
regularize 0 O in B using the equation described in Remark 
3. This type of approach is addressed in the work of Chan 
and Shen mentioned previously [9], 

5) In practice, functional (10) is used to interpolate shapes, 
i.e., to interpolate level sets. The image is decomposed 
into upper level sets [u 0 > A], which are interpolated 
using (10) to produce the level sets X\u of a function 14, 
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which is reconstructed inside 0 by using the reconstruc- 
tion formula 

u(x) = sup{A : x e X\u}. 

To guarantee that the reconstructed level sets correspond 
to the level sets of a function u, they should satisfy that 
X\+iu C X\u. In practice, we force our solution to 
satisfy this property. 

In principle, our functional (10) could be used directly 
to interpolate functions. But, discontinuities of the image 
have a contribution to the energy which is proportional 
to the jump. This gives different weights to discontinu- 
ities of different sizes and, as a consequence, they are 
not treated in the same manner. This is not reasonable if 
we want to interpolate the shapes of the image, indepen- 
dently of their contrast. When taking level sets, we treat 
all shapes equally, and the parameters of the functional 
should only weight geometric quantities (like length, total 
curvature) and decide which interpolation is taken as a 
function of them. This approach is less diffusive than di- 
rectly interpolating the gray levels. Theorem 2 in Ap- 
pendix I proves the existence of minimizers for our model 
and can be applied to both cases, binary and gray level 
images. It guarantees that there are minima of our func- 
tional. We do not yet know the qualitative properties of 
those minima. Even if at the intuitive level the main quan- 
tities that appear in our functional are length and a mea- 
sure of total curvature (like in the elastica), we do not have 
a rigorous proof that this is so. The functional was in- 
troduced on a heuristic basis, but relaxations may occur 
as they occur in (16), where div(Vu/|Vu|) may repre- 
sent perimeter/area when computed on a flat region [3], 
[37]. This requires further study and we shall pursue it 
elsewhere. 

6) The choice made in Remark 5 of decomposing the image 
uq into upper level sets, interpolating them and recon- 
structing the function u, introduces a lack of symmetry. 
Indeed, we are giving more weight to upper level sets than 
to lower level sets. This can be seen in Fig. 3. Fig. 3(a) 
displays the image to be interpolated. It is clear that sev- 
eral reasonable solutions are possible and no one of them 
is preferable to the others. The choice we made gives 
Fig. 3(b) as solution, favoring that the object whose level 
is 2 1 0 goes above the object whose level is 0. But, in that 
case, the "true" information is lacking and we selected 
one of the possible reasonable solution. 



C. Interpolating Gray Values along the Integral Curves of a 
Vector Field 

Our purpose in this section is to further discuss the term 

\Vu\-0-Vu (17) 



L 





1 

% £ . ■ I i ill! 



(a) 



(b) 



Fig. 3. (a) Original image and (b) result of fUling-in with our algorithm and 
the level-sets ordering selected in this paper. 

from image gradients). We shall see that (when 0 is known), 
when minimizing (17), we are constructing the function u 
whose values on the boundary are given and whose direction of 
the gradient is given by 0. We shall discuss this from a general 
point of view. Thus, suppose that Q, is an open bounded domain 
with a Lipschitz boundary and \p e L°°(dtt). Let // : Q -> Er 
be a vector field whose smoothness will be detailed below. We 
ask the following question: can we interpolate the boundary 
data ip along the integral curves of vl In the case discussed in 
last section, v = 0- 1 , and we propagate the boundary data tp 
along the integral curves of 0 X . Heuristically, 0 is orthogonal to 
the level lines of u. Coming back to our general discussion, we 
want to construct a function u : Q —> M such that u\oq — <p 
and u being constant along the integral curves of //, i.e., the 
solutions of the system of ordinary differential equations 



dX 



This amounts to say that 



v • Vu = 0 



(18) 



(19) 



in functional (10) (see also [24] for a related, L 2 and 
Poisson-equation based, approach of gray value reconstruction 



a first-order transport equation whose characteristic curves are 
the solution of the system (18). Let us discuss the difficulties 
posed by this formulation. First of all, existence and uniqueness 
of solutions of (18) is guaranteed when v is a Lipschitz vector 
field, a very strong regularity assumption, which excludes any 
singularity for v. More general existence results have been ob- 
tained in [13] via the study of transport equations, indeed, via 
formulations analog to (19). Typically, they are assuming that // 
is in some Sobolev space like e (M N ), with some other in- 
tegrability assumptions, and div(z/) 6 L°°(JR N ). These results 
have been further extended in [12], [30]. In particular, Lions in 
[30] proves a.e. existence of solutions of (18) for vector fields 
which are piecewise in W 1,1 in a precise sense defined by the 
author. As observed in these papers, it is not known if the pre- 
vious result is true for BV vector fields. On the other hand, 
even for a vector field in W 1 * 1 , for which we have existence 
a.e. of solutions of (18), the problem of constructing u satis- 
fying (19) and such that u\oq = cp is not obvious. Indeed, con- 
sider a smooth vector field v defined on a simple domain, like 
D := {x £ M 2 : \x\ < 1) and suppose that the integral curves 
of u are curves that foliate D and such that at any point of 3D 
we start a curve that ends in another point of 3D. Then the only 
possibility to extend u to fl so that u\on = in a classical sense 
is that <p takes the same values at the beginning and endpoints 
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Fig. 4. (a) Original image, (b) original image with missing region, and (c) Filled-in image when the gradients inside the white square are known. 

of the integral curves of v A possibility to overcome this dif- 
ficulty, would be to use the vanishing viscosity method, i.e., to 
solve the elliptic equation 

v ■ Vu e + eAu € = 0 (20) 

and let e — ► 0. Then, we hope the sequence u € to converge to 
some bounded function u which solves the problem in a distri- 
butional sense. We do not have further information on the reg- 
ularity of u. On the other hand we do not know in which sense 
the boundary conditions hold. 

Let us consider the problem from the algorithmic point of 
view, i.e., we want to design an effective algorithm to solve it. 
Since the problem may be ill-posed, because of incompatibility 
of boundary data joining two integral curves of the vector field 

we propose a variational formulation of the problem. Let 0 = 
u 1 - . Assume that |0| < 1. If a solution exists, then 0 should 
point in the normal direction to the level lines of u. We implicitly 
assume that 9 should be constructed as the vector field normal 



to the level curves of u. Then, formally, 6 • S7u — |Vu|. Thus, 
it seems reasonable to minimize the functional 



I'Xu) = [ |Vu| - /V Vu 
Jq Jq 



(exactly the one introduced above) defined in the set of functions 
of bounded variation BV(£2) whose trace at the boundary is 
given by <p. Let us formally integrate by parts in the second term 
of F(u) to obtain 

F(u)= [ |V-u|+ f div(0) u- / 0'™' 
Jq Jq Jon 

Since u y 0 are known at the boundary, minimizing i r amounts 
to minimize 

E(u) = / |V«| + / div(0) • u. 

Let us make precise the class of admissible functions where E is 
minimized. We assume that div(0) € and <f € L°°(8Q). 
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It seems reasonable to impose that the solution u is a bounded 
function with an L°° bound given by |M|<x> (or a constant re- 
lated to IMloo and the size of ft). Then the second integral in the 
definition of E(u) is well defined. The first integral requires the 
use of the space of bounded variation functions. Thus our admis- 
sible class is A = {u € BV(ft) : \u(x)\ < \\y>\\oo a.e. u\ou = 
(/?}. We propose 



Minimize / \Vu\ + / div(#) • u 
Jn Jn 



u e A. 



(21) 



As is well-known [1 5], [21], the solution of this problem has to 
be understood in a weak sense as the solution of the problem 



Minimize / \Vu\ 4- / div(#) - u+ i \u 

Jq Jq Jon 

u 6 BV(ft) 

M < li<p||oo. 



/ |« - <p\ dH l 
Jon 



(22) 



Then we have the following result. 

Theorem 1: Let 0 e L} oc (ty 2 , with div(0) e L 1 ^),^ € 
L°°(dft). Then there is a function u e BV(ft) such that 
\ u ( x )\ 5: IMIoc ae - minimizing (22). 

Proof: The result is contained in [2 1], Theorem 1 .4. □ 

This clarifies the role of the term (17) in (10). 

III. Numerical Experiments 

To minimize (21) we use the steepest descent method. For 
that, we formally compute the Euler-Lagrange equation for u, 
namely, 

~ div (j|^)+ div W = 0 ( 23 > 

supplemented with Dirichlet boundary conditions for u. In prac- 
tice, we use the evolution equation 



u t 



= div ^ 



\Wu\) 



div(O) 



with Dirichlet boundary data and initial condition constructed as 
an ad-hoc interpolation that will be corrected by the equation. 
General existence results which can be adapted to this case can 
be found in [4]. Note that the vector field 0 is assumed to be 
known in this case. This limits the usefulness of this model. But 
we present some experiments below to illustrate the role of this 
term. 

To minimize the functional (10) we use the steepest descent 
method. If we denote the energy term by E(9. u), the steepest 
descent equations are 



and 



0 t = -V e E{9,u) in ft 
Ut = -V u E(9,u) in a 



(24) 
(25) 



supplemented with the corresponding boundary data and initial 
conditions. The constraints on (9, u) can be incorporated either 
by penalization or by brute force after each time step. Given 
X C JR 2 we denote by xx the characteristic function of X, 



J £ 

^jgj mm? 




(a) (b) 
Fig. 5. Simple example of our algorithm. 





(a) (b) 
Fig. 6. Demonstration of the freedom in topology of our algorithm . 

i.e., Xx{x) = 1 if ^ € X, otherwise, xx{%) = 0- To simplify 
our notation, let us write g = b\div(9)\ p , h — a -f b\Vk * u\. 
Then 

V e E(9,u) = -pV[h\&v(0)\ p - 2 div(0)] 

- a(Vuxn + VuoXb) = 0 (26) 

and 

V. £ <»„.) = - d ,v(,.(^))-„ d ,v(^) 

+ adiv(0) = O. (27) 

In our experiments, we take k a Gaussian kernel with small vari- 
ance, say one or two pixels. In practice, one can also dismiss 
the kernel k. These equations have to be complemented with 
the corresponding boundary conditions for 9 and u specified by 
the admissible class, i.e., we specify the normal component of 
0 in Oil and the Dirichlet boundary condition for u in di 1, since 
u = uo in B. The initial conditions are ad-hoc interpolations, 
for instance, we can take u inside ft as the average value of Uq 
in B, 9 inside ft being the direction of the gradient of u. One 
can also take a geodesic propagation inside ft of the values of 
uo in B, with 9 being again the direction of the gradient of u. 
The exact numerical implementation of these flows is given in 
Appendix II. 

In the experiments below, this algorithm is used to interpolate 
level sets, following the approach in [31] and [32]. The image 
in B is decomposed into level sets and we get a family of binary 
images u$x = X[a Q >x]A — 0, 1,2, ... , 255. These functions 
are interpolated inside ft and we obtain a family of level sets 
X\u. Then the function u is reconstructed using the reconstruc- 
tion formula 

u(x) = sup{A e {0, 1, ... , 255} : x € X x u}. 

As observed in Remark 5 of Section II-B, we force our solution 
to satisfy the monotonicity property of the level sets, i.e., that 
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Fig. 7. Example of the proposed filling-in algorithm. 




Fig. 8. Detail of the example in Fig. 7. Note the smooth continuation of the 
edges. 




(b) (c) 

Fig. 9. Level-set corresponding to the region in Fig. 8. (a) Original and the 
results with (b) p = 1 and (c) p — 2 . 

X\+\u C Xya. This is imposed in the initialization of the level 
set X\u and is maintained at each iteration of the algorithm by 
taking the supremum of the current solution with the charac- 
teristic function of X\u. With this approach, we diminish the 
diffusive effects of the above algorithm and we better capture 
the shapes and discontinuities on the interpolated image. 

The constraints on 0 and ||u||oo can be introduced after each 
iteration of the above equations. We also comment that the con- 
straint 0 - Vu = |V'ti|, which was introduced as a penalization 
term, could also be introduced by brute force after each time 
step iteration of the algorithm. Let us describe the experiments. 




Fig. 10. Example of automatic text removal. 



First, in Fig. 4 we display some experiment to illustrate func- 
tional (21). Fig. 4(a) displays the full image without the hole. 
Fig. 4(b) displays the image with the hole. The vector field 9 has 
been computed on Fig. 4(a) and we see in Fig. 4(c) the result of 
interpolating the gray level knowing the vector field inside Q. 
We see that the shape of the eye is recovered but not the gray 
level. This is not a surprise since the gray level inside the eye 
cannot be recovered from the gray level on the boundary of Q. 
The algorithm is able to capture the shapes inside the eye by in- 
tegrating the vector field 0. 

In the following experiments, we show the results of the joint 
interpolation of gray level and the vector field of directions using 
functional (10). The experiments have been done with p = 1 
and/or p = 2. The results are quite similar. Unless explicitly 
stated, we display the results obtained with p = 1. Fig. 5(a) 
displays an image made of four circles covered by a square. In 
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Fig. 1 1. Example of automatic graffiti removal 



Fig. 5(b) we display the result of the interpolation. In Fig. 6(a), 
we display an example where the hole is not simply connected. 
The interpolation is displayed in Fig. 6(b). Fig. 7(a) is the image 
of Lena with two holes, a lower one in the hat and an upper one. 
Fig. 7(b) displays the result of the interpolation. Fig. 8 displays a 
zoom of the region around the lower hole. In Fig. 9(a) we display 
a level set of uq corresponding to the region around the lower 
hole. Fig. 9(b) and (c) display the corresponding interpolation 
withp = 1 andp = 2, respectively. Fig. 10(a) displays an image 
with text to be removed. Fig. 10(b) displays the corresponding 
reconstruction result. Fig. 1 1(a) displays a portion of an image 
with text. Fig. 1 1(b) displays the corresponding reconstruction 
result, obtained with p = 2. 

IV. Concluding Remarks 

In this paper, we have proposed a formal variational approach 
for filling-in regions of missing data in still images. The basic 
idea is to smoothly extend inside the hole both the vector field 
obtained from the image gradient and the corresponding gray 
values. We have presented a number of examples and showed 
theoretical results regarding the proposed formulation. 

A number of research directions are suggested by the work 
here presented. First of all, we need to complement this algo- 
rithm by a technique capable of filling-in textured regions. Sec- 
ondly, the extension of the framework to the filling-in of other 
type of missing imagery data is of great interest for a number of 
applications. Last, we would like to study these ideas for inter- 
polation in video data. These topics will be the subject of sub- 
sequent reports. 

Appendix I 
Existence of Minimizers 

Recall that ft = ft U B is an open bounded set 
whose boundary is Lipschitz. For simplicity, let us de- 
Fine the class B of admissible pairs (0,'u) where 0 e 

W l *(div,n),u £ BV(ft,£,uo),|0| < M«| < il'"o|L~<a) 
md9n\ on = 0 o ■ n| 0 „. 

Theorem 2: If p > 1, there is a minimum (0, u) e 3 for the 
problem (10). 

Proof: Let us denote by E(6,u) the energy defined in 
(10). Let (0 nj « n ) be a minimizing sequence for E(0, u). Since 



L 



\vu n \-e. n - vix n >o 



and E(9 n .u n ) is bounded, we obtain that 

idiv(0„)r 



is bounded. Since 1 9 n \ < 1, we have that 9 n is weakly relatively 
compact in all spaces L q (ft) 2 for all 1 < q < oo and we may 
assume that 0 n — > 9 weakly in L q (ft) 2 for all 1 < q < oo and in 
PF 1)P (div, ft). Now, integrating by parts the term fa 9 n - Vu n , 
we obtain 

/ 0 n • Vu n = " / div(0 n )tt„ + / [9nMUn 

JCi Jq Jon 

= - / div(9 n )u n -h f [0o,n]u o . 

JU JOQ 

The integration by parts is possible by results of Anzellotti [2] 
given above. From the above identity, we obtain 

l/Vvun <||div(fl n )IWKII^+ / Kl 

\Jq Jon 

since u n = u 0 in dft, where p' is the exponent conjugated to 
p. Since we minimize the energy E(9,u) for functions with an 
L°° bound, we obtain that 



9 n -VUn 



is uniformly bounded in n. The consequence of this observation 
is that 



I*} 



|V« n | 



L 



is also bounded. Then, modulo a subsequence, we may assume 
that u n converges to some function u in L l (ft). Note that u € 
B V( ft, B, uq). Since we have an L°° bound on u Tl , we also have 
that u n converges to u in L q (ft) for all 1 < q < oo. Then 
Vfc * u n — > V/c * u uniformly in ft. In particular, we obtain 

/ |div(e)| p (a + 6|V**ti|)£te 
Jn 

<liminf / \di\(O n )\ p {a + 6| Vk * u n \) dx 
" Jn 

and 

/ IVlil < lim inf / IVitnl. 
Jn 11 Jn 

Finally, since div(0 T ,) weakly converges to div(^) in LP (ft) and 
u n — > u in LP'(ft), passing to the limit in 

/ 9 n • Vit^ = - / div(6 n )u n + / [8o 7 n]uo 
Jn Jn Jon 

we get that fa 6 n • Vu n converges to 

- / div(0)u + / [0 0 ,n]u o - fo-X/u. 
Jh Jon Jn 
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Thus, collecting all these facts, we obtain that 
E{0,u) < Km'mfE(0 n) u n ). 

The pair (6>u) is a minimum of E in the class of admissible 
functions for this functional. □ 

Appendix II 
Numerical Implementation 

To solve (24) and (25), we use a implicit discretization in 
time. To be precise, we write 

V*E(M',t* : v) - -pV[A(e+ idiv(^|*- 2 )div(0)] 

-a{Vu X n + VuoXs) =0 (28) 

and 

_ ' >< " v (^TCT) +<,div( '' , 

= 0. (29) 
Then, discretize in time according to 

0 n + l -6 n = &tV e E{e n +\0 n ,u n ,u n ), (30) 

and 

u »+i _ u n = AtV tt f?(» u+1 ,» w+l J ti n+1 ,tt")- ( 31 ) 
Finally, we make the change of variables £ n+1 = $ n+1 - 

= A*V # £(r +1 + 0 n , (32) 
= A*V u JE(^ + S0 7l+ \v" +l +u",it")- ( 33 ) 

Now, since 6 n -n^ = f? n+1 -n)^ and ?/ n |an = u n+l |0O. then 
the normal component of f n+1 and the value of v n+l are zero 
on the boundary, and we may use a conjugate gradient method 
to solve (32) and (33). The constraint |0| < 1 is incorporated 
by brute force after each time step. We can also set a = 0 and 
incorporate the constraint that |Vw| — 6 ■ Vu by brute force 
after each time step. 
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